Load {sf} and read in a subset of the NH National Hydrography Dataset (available from www.granit.unh.edu):

#install.packages("sf")
library(sf)

local_lakes <- st_read("./data/NHDWaterbodyDataset_NHSeacoast.shp")

Load {sf} and read in a subset of the NH National Hydrography Dataset (available from www.granit.unh.edu):

#install.packages("sf")
library(sf)

local_lakes <- st_read("./data/NHDWaterbodyDataset_NHSeacoast.shp")
class(local_lakes)

Load {sf} and read in a subset of the NH National Hydrography Dataset (available from www.granit.unh.edu):

class(local_lakes)

## [1] "sf"         "data.frame"

Load {sf} and read in a subset of the NH National Hydrography Dataset (available from www.granit.unh.edu):

head(local_lakes)

Load {sf} and read in a subset of the NH National Hydrography Dataset (available from www.granit.unh.edu):

head(local_lakes,3)

## Simple feature collection with 3 features and 13 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -71.28925 ymin: 42.91414 xmax: -70.7096 ymax: 43.44024
## Geodetic CRS:  NAD83
##     Prmnn_I      FDate Resoltn GNIS_ID GNIS_Nm    AreSqKm Elevatn
## 1 141035025 2005-10-04       2    <NA>    <NA> 0.00100000      NA
## 2 141036768 2019-03-19       2    <NA>    <NA> 0.00028992      NA
## 3 141032014 2019-03-19       2    <NA>    <NA> 0.00038662      NA
##          ReachCd FType FCode VsbltyF      Shp_Lng      Shap_Ar
## 1 01060003011689   390 39004   24000 0.0010302955 7.489623e-08
## 2           <NA>   466 46600       0 0.0007273591 3.203682e-08
## 3 01060003008954   390 39004   50000 0.0007978068 4.298463e-08
##                         geometry
## 1 POLYGON ((-70.86617 42.9143...
## 2 POLYGON ((-71.28921 43.0608...
## 3 POLYGON ((-70.70971 43.4399...

Interactive mapping tools are great for sharing and communicating spatial data, but also when you just want to look at the data.

One of these packages is {mapview}. Let’s take a look at our local_lakes dataset:

#install.packages("mapview")
library(mapview)

mapview(local_lakes)

{sf} integrates well with ‘tidy’ data workflows. What do I mean by that?

{sf} integrates well with ‘tidy’ data workflows. What do I mean by that?

Let’s subset our local_lakes dataset to only include my field site, Swains Lake:

library(dplyr)

swains_lake <- local_lakes %>%
  filter(GNIS_Nm == "Swains Lake")

{sf} integrates well with ‘tidy’ data workflows. What do I mean by that?

Let’s subset our local_lakes dataset to only include my field site, Swains Lake:

library(dplyr)

swains_lake <- local_lakes %>%
  filter(GNIS_Nm == "Swains Lake")

mapview(swains_lake)

The st_read() function automatically stores information about our data, including:

• the data format
  
• the coordinate reference system (CRS)
  
• the spatial extent of the data

The st_read() function automatically stores information about our data, including:

• the data format
• the coordinate reference system (CRS)
  
• the spatial extent of the data

st_geometry_type(swains_lake)

## [1] POLYGON
## 18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT ... TRIANGLE

The st_read() function automatically stores information about our data, including:

• the data format
  
• the coordinate reference system (CRS)
  
• the spatial extent of the data

st_crs(swains_lake)

## Coordinate Reference System:
##   User input: NAD83 
##   wkt:
## GEOGCRS["NAD83",
##     DATUM["North American Datum 1983",
##         ELLIPSOID["GRS 1980",6378137,298.257222101,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["latitude",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["longitude",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4269]]

Proper representation of coordinate reference systems (CRS) is very important for all spatial work!

CRS are often denoted using an “EPSG” code (see bottom of st_crs output). Common CRS include:

• WGS84 (EPSG = 4326)
  
• NAD83 (EPSG = 4269)

See other resources for more information, including https://datacarpentry.org/organization-geospatial/03-crs/

st_crs is your friend.

st_crs(swains_lake)

The st_read() function automatically stores information about our data, including:

• the data format
  
• the coordinate reference system (CRS)
  
• the spatial extent of the data

st_bbox(swains_lake)

##      xmin      ymin      xmax      ymax 
## -71.05927  43.18261 -71.02456  43.19972

The st_read() function automatically stores information about our data, including:

swains_lake

## Simple feature collection with 1 feature and 13 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -71.05927 ymin: 43.18261 xmax: -71.02456 ymax: 43.19972
## Geodetic CRS:  NAD83
##     Prmnn_I      FDate Resoltn  GNIS_ID     GNIS_Nm  AreSqKm Elevatn
## 1 141033653 2015-01-14       2 00870298 Swains Lake 1.379346      NA
##          ReachCd FType FCode VsbltyF   Shp_Lng      Shap_Ar
## 1 01060003002342   436 43600 1000000 0.1605085 0.0001527365
##                         geometry
## 1 POLYGON ((-71.05117 43.1887...

Say we want to find the closest weather station to Swains Lake.

Say we want to find the closest weather station to Swains Lake.

We can use the {rnoaa} package to identify NOAA climate stations, and then use spatial joins to select the nearest site and download some data.

library(rnoaa)

# Takes a minute to run:
stations <- rnoaa::ghcnd_stations() %>% 
  filter(state=="NH",element=="PRCP",last_year>2020) %>%
  st_as_sf(.,coords=c("longitude","latitude"),crs=4326) %>%
  arrange(id)

stations <- read_sf("./data/NH_GHCND_stations.shp") 

Say we want to find the closest weather station to Swains Lake.

We can use the {rnoaa} package to identify NOAA climate stations, and then use spatial joins to select the nearest site and download some data.

head(stations,3)

## Simple feature collection with 3 features and 9 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -71.55 ymin: 43.4356 xmax: -71.3231 ymax: 43.5541
## Geodetic CRS:  WGS 84
## # A tibble: 3 × 10
##   id          elevation state name  gsn_flag wmo_id element first_year last_year
##   <chr>           <dbl> <chr> <chr> <chr>    <chr>  <chr>        <int>     <int>
## 1 US1NHBK0001      207. NH    TILT… <NA>     <NA>   PRCP          2009      2021
## 2 US1NHBK0002      214  NH    BELM… <NA>     <NA>   PRCP          2009      2021
## 3 US1NHBK0007      172. NH    LACO… <NA>     <NA>   PRCP          2009      2021
## # … with 1 more variable: geometry <POINT [°]>

Calculate the distances between each weather station and Swains Lake.

st_distance(stations,swains_lake)

Calculate the distances between each weather station and Swains Lake.

We need to transform the stations and swains_lake data so that they share the same CRS.

# Transform data to CONUS Albers Equal Area Conic projection:
swains_lake_proj <- st_transform(swains_lake,5070)
stations_proj <- st_transform(stations, 5070)

Calculate the distances between each weather station and Swains Lake.

stations_proj2 <- stations_proj %>%
  mutate(dist = st_distance(.,swains_lake_proj,by_element = TRUE))

stations_proj2 %>% arrange(dist) %>% select(id,name,dist) %>% head(3)

Spatial join: Find the closest weather station to Swains Lake.

closest_station <- st_join(swains_lake_proj,stations_proj,join=st_nearest_feature)

Spatial join: Find the closest weather station to Swains Lake.

closest_station <- st_join(swains_lake_proj,stations_proj,join=st_nearest_feature)
head(closest_station)

## Simple feature collection with 1 feature and 22 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 1993730 ymin: 2506799 xmax: 1996588 ymax: 2508389
## Projected CRS: NAD83 / Conus Albers
##     Prmnn_I      FDate Resoltn  GNIS_ID     GNIS_Nm  AreSqKm Elevatn
## 1 141033653 2015-01-14       2 00870298 Swains Lake 1.379346      NA
##          ReachCd FType FCode VsbltyF   Shp_Lng      Shap_Ar          id
## 1 01060003002342   436 43600 1000000 0.1605085 0.0001527365 US1NHST0040
##   elevation state             name gsn_flag wmo_id element first_year last_year
## 1      80.8    NH BARRINGTON 3.2 E     <NA>   <NA>    PRCP       2016      2021
##                         geometry
## 1 POLYGON ((1994469 2507096, ...

Now we can download and plot our precipitation data:

prcp_data <- rnoaa::ghcnd_search(closest_station$id, var = "PRCP")

# Plot data:
prcp_data$prcp %>% filter(date>"2021-08-01" & date<"2021-10-01") %>% 
  mutate(precip_mm = prcp/10) %>%
  ggplot() + geom_line(aes(x=date,y=precip_mm)) +
  labs(x = "Date (2021)", y = "Precipitation (mm)") + theme_classic()

Useful packages for mapping in R:

• {ggplot2} along with geom_sf()
  
• {tmap} for customized, publication-quality maps
  
• {mapview} and {leaflet} for interactive maps!
  
• {rasterVis} for raster data

We can still use {ggplot2} for plotting spatial data.

ggplot() + 
  geom_sf(data=swains_lake) + 
  coord_sf()

We can still use {ggplot2} for plotting spatial data.

A lot of the formatting should already be familiar if you’re a {ggplot2} user:

ggplot() + 
  geom_sf(data=swains_lake,fill="cyan3",color="cyan4") + 
  coord_sf(ylim = c(43.18, 43.2)) + 
  ggtitle("Swains Lake, NH") + 
  theme_classic() 

Let’s load one more spatial dataset to include in our map - a shapefile of the watershed boundary, that is, the area of land that drains into Swains Lake.

swains_watershed <- read_sf("./data/Swains_watershed_boundary.shp")
st_crs(swains_watershed)

## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["latitude",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["longitude",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4326]]

Let’s load one more spatial dataset to include in our map - a shapefile of the watershed boundary, that is, the area of land that drains into Swains Lake.

ggplot() + 
  geom_sf(data=swains_lake,fill="cyan3",color="cyan4") + 
  geom_sf(data=swains_watershed, fill=NA, color="gray20", lwd=0.6) +
  coord_sf() + 
  ggtitle("Swains Lake, NH") + 
  theme_classic() 

Add scales, orientation arrows, etc. using {ggspatial}.

You can probably think of other ways you want to edit your map for a report or publication, but we’re definitely getting closer!

library(ggspatial)

swains_map <- ggplot() + 
  annotation_map_tile(zoom = 15) +
  geom_sf(data=swains_lake,fill="cyan3",color="cyan4") + 
  geom_sf(data=swains_watershed, fill=NA, color="gray20", lwd=0.6) +
  ggtitle("Swains Lake, NH") + 
  theme_classic() +
  labs(x="Longitude",y="Latitude") + 
  coord_sf() +
  annotation_north_arrow(location="bl", 
                         height = unit(.75, "cm"),
                         width = unit(.5, "cm")) +
  annotation_scale(location="bl", 
                   pad_x = unit(1, "cm"), 
                   text_cex = 0.5, 
                   height = unit(0.2, "cm")) 

Add scales, orientation arrows, etc. using {ggspatial}.

You can probably think of other ways you want to edit your map for a report or publication, but we’re definitely getting closer!